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Abstract 


The subfamily Plectorhinchinae (sweetlips) is composed of poorly-known species with high commercially and ecologically val- 
ues that exhibit phenotypic plasticity and various morphologies. Few studies have assessed the validity of sweetlips, intergeneric 
relationships and evolutionary survey in this subfamily, which have not yet been resolved. This study investigated the DNA se- 
quences of (1) the mitochondrial COI gene to delimit species, and (2) two mitochondrial (COI and Cyt b), and one nuclear (RAG1) 
markers to infer phylogenetic relationships and evolutionary and biogeographic history. The molecular results could differentiate 
Diagramma punctatum from the other species, but failed to distinguish D. /abiosum as a distinct species with considerably lower 
genetic distances for the COI (0.53%) and Cyt b (0.51%) markers. However, additional taxonomic investigations are required to 
shed light on this issue. All previously described nominal species of sweetlips in the northwest Indian Ocean were found to be well 
supported. The monophyly of Plectorhinchus is not supported and Diagramma pictum and D. punctatum should be assigned to the 
genus Plectorhinchus. The biogeographic history of Plectorhinchinae likely originated in the Indo-Pacific ca. 34 Ma (30-39 Ma; 
late Eocene/ middle Oligocene) and subsequently colonised the Western Indian Ocean and the Central Indo-Pacific. Maximum di- 
versification within the subfamily occurred from the middle Miocene to Pliocene, coinciding with dispersal and vicariance events. 
Diversification was probably driven by both biological and geographical factors. 
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Introduction 21 genera and 136 species (Fricke et al. 2022). The Plec- 
torhinchinae commonly known as sweetlips, inhabit the 


Accurate species delimitation and phylogenetic recon- — shallow coastal waters of the Indo-West Pacific and East- 


struction are vital to understand biodiversity assessments, 
conservation management, evolutionary patterns, and 
processes (Traldi et al. 2020; McCord et al. 2021). Mem- 
bers of the family Haemulidae Gill, 1885 are dominant in 
shallow coral or rocky reefs of tropical and subtropical 
waters worldwide (McKay 2001; Lindeman and Toxey 
2003; Nelson et al. 2016). They also have commercial and 
ecological values for most nations (Rocha et al. 2008). 
The family comprises two subfamilies Haemulinae Gill, 
1885 and Plectorhinchinae Jordan & Thompson, 1912, 


ern Atlantic, which currently comprises 42 recognised 
species (Fricke et al. 2022). The species are classified into 
four genera (Diagramma Oken, 1817, Genyatremus Gill, 
1862, Parapristipoma Bleeker, 1873 and Plectorhinchus 
Lacépede, 1801), of which Plectorhinchus is the most 
speciose (Fricke et al. 2022); however, previous molecu- 
lar studies have shown that Genyatremus does not belong 
to the Plectorhinchinae (Bernardi et al. 2008; Sanciangco 
et al. 2011; Tavera et al. 2018). Ever since the validation 
of Diagramma and Plectorhinchus, the identification of 
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haemulid species, including several cases of re-descrip- 
tions and re-allocations has been intensively debated 
(Johnson et al. 2001; Johnson and Wilmer 2015; Tavera et 
al. 2018; Damadi et al. 2020). In reality, further research 
on phylogenetic relationships and taxonomic treatment 
which, in the past, has been less identified in these genera 
of Haemulidae is required to figure out the actual num- 
ber of species because of nominal and complex species 
as well as the intergeneric relationships between these 
groups of fishes (Johnson and Wilmer 2015; Tavera et 
al. 2018). Moreover, many haemulid species manifest- 
ed morphological diversity in colouration and markings 
with developmental stages, for example, Plectorhinchus 
species, or population variation associated with changes 
in environmental conditions, such as Diagramma pictum 
(McKay 1983; Liang et al. 2016). Thus, adopting a mul- 
tiple analytical approach to DNA sequence data can be 
promising (Padial et al. 2010; Carstens et al. 2013) be- 
cause it evaluates the morpho-molecular diversity to un- 
derstand their phylogenetic relationships and lowers the 
possibility of homoplasy (Nei and Kumar 2000; Tavera 
et al. 2011; Varon-Gonzalez et al. 2020). The genera 
Diagramma and Plectorhinchus comprise 35 valid spe- 
cies worldwide (Fricke et al. 2022). Recent molecular 
studies confirmed that the subfamily Plectorhinchi- 
nae 1s monophyletic, but the phylogenetic relationships 
amongst some genera are challenged by those retrieved 
by morphological studies (Sanciangco et al. 2011; Tavera 
et al. 2018). Johnson et al. (2001) examined the genus 
Diagramma and five geographically constrained subspe- 
cies recognised within D. pictum (Thunberg, 1792), based 
ona morphological approach: D. p. pictum from Indo-Ma- 
lay Archipelago to Central Indo-Pacific excluding Austra- 
lia (type locality: Japan); D. p. labiosum from northern 
Australia to South New Guinea (type locality: Australia, 
Queensland); D. p. cineraseens from the Persian Gulf to 
the Bay of Bengal (type locality: Sri Lanka, Trincomalee); 
D. p. punctatum from the Red Sea (type locality: north- 
ern Red Sea) and D. p. centurio occurring along the 
African coast (type locality: Seychelles). The taxonom- 
ic position of Digramma members has been challenged 
(Allen 1997; Johnson et al. 2001; Heemstra and Heemstra 
2004; Bogorodsky et al. 2014; Johnson and Wilmer 2015; 
Bray 2017; Golani and Fricke 2018). Recently, many 
subspecies have been promoted to the status of species 
rank (Fricke et al. 2022). After examining several spe- 
cies, Tavera et al. (2018) proposed to classify them under 
Plectorhinchus, but it is still necessary to investigate the 
remaining nominal species (D. labiosum Macleay, 1883, 
D. punctatum Cuvier, 1830 and D. melanacrum Johnson 
& Randall, 2001). Of 42 species, 11 species are morpho- 
logically known from the northwest Indian Ocean, includ- 
ing Diagramma pictum, D. punctatum, Plectorhinchus 
albovittatus (Ruppell 1838), P flavomaculatus (Cuvier, 
1830), P playfairi (Pellegrin, 1914), P gaterinus 
(Forsskal, 1775), P. gibbosus (Lacepéde, 1802), P. pictus 
(Tortonese, 1936), P. sordidus (Klunzinger, 1870), 
P. schotaf (Forsskal, 1775) and P. makranensis Damadi, 
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Yazdani Moghaddam, Ghassemzadeh & Ghanbarifar- 
di, 2020. A comprehensive analysis of all species in this 
region has not been done. Some of these species have 
limited distribution, such as D. punctatum (Red Sea) and 
Plectorhinchus makranensis (the Persian Gulf to the Gulf 
of Oman) (Zajonz et al. 2019; Damadi et al. 2020). In 
addition, few genetic studies have explored haemulid 
species, including DNA barcoding analysis in a few taxa 
from the northwest Indian Ocean (Asgharian et al. 2011; 
Rabaoui et al. 2019). The main goals of this study are as 
follows: (1) to test the validity of the northwest Indian 
Ocean taxa proposed by morphological classification, (2) 
to utilise the wide geographic coverage of sweetlips spe- 
cles using species-delimitation algorithms to shed light on 
their molecular operational taxonomic units (MOTUs), 
(3) to reconstruct robust phylogenetic trees using two 
mtDNA (COI and Cyt 6 and one ncDNA (RAGI1) markers 
to infer the phylogenetic relationships and the systematic 
position of Diagramma pictum and D. punctatum and (4) 
to estimate divergence times between genera and species 
and reconstruct ancestral ranges to test the relative impor- 
tance of historical geology, biology and ecological factors 
for the diversification within major clades of the subfam- 
ily Plectorhinchinae. 


Materials and methods 
Sample collection and specimen examination 


A total of 254 individuals belonging to 29 species, in- 
cluding both our new samples and the sequences taken 
from GenBank and the Barcode of Life Data systems 
(BOLD systems) were analysed (Suppl. material 1). A to- 
tal of 44 new sequences of three molecular markers were 
used. The specimens were photographed and then a small 
piece of right pectoral fin of each individual was excised 
and preserved in absolute alcohol. The whole samples 
were stored in formaldehyde for long-term storage. All 
collected fish specimens were identified as described by 
Randall (1995). Voucher specimens are deposited in the 
Zoological Museum, Ferdowsi University of Mashhad 
(ZMFUM), Iran (Suppl. material 1). 


DNA extraction, amplification and sequencing 


Genomic DNA was extracted from the muscle tissue 
following the BioGene kit protocol. The reactions volume 
of polymerase chain reaction (PCR) was 25 ul, including 2 
ul template DNA, 9.5 yl ddH,O, 12.5 pl Master Mix and 
0.5 ul of each forward and reversed primer with 10 uM 
concentration. Two mitochondrial and one nuclear gene 
were amplified as follows: 617 bp of the CO1 was ampli- 
fied using universal primers COILBC_F and COIHBC_R 
(Ward et al. 2005), 1100 bp of the Cyt 6 was amplified using 
UniF and UniR primers (Orrell et al. 2002) and 1440 bp of 
the nuclear gene RAG1 was amplified using nested PCRs 
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product (1** PCR: 2510F and RAGIR1; 2"? PCR: RAGIF1 
and RAGIR2 primers) (Lopez et al. 2004). PCR conditions 
for the COI, Cyt 6, and Ragl genes were as reported by 
Sanciangco et al. (2011). The PCR products were evaluat- 
ed with 1.5% agarose gel for band characterisation, puri- 
fied and then sent for sequencing at Microsynth Company, 
Switzerland. Ocyurus chrysurus (Bloch, 1791) and Rhom- 
boplites aurorubens (Cuvier, 1829) from the family Lut- 
janidae Gill, 1861 were used as outgroups, since Lutjanid 
species has a close phylogenetic relationship with Haemul- 
idae (Betancur-Rghton et al. 2013; Near et al. 2013). 


Phylogenetic analysis and population structure 


The quality of new sequences was probed using the 
CLC Genomics Workbench v.3.6.5 software (QIAGEN 
Bioinformatics). The sequences were aligned in MAFFT 
v.7.463 and edited with BIOEDIT v.7.2 software (Hall 
1999). Genetic interspecific distances of the mitochondrial 
COI gene were calculated by MEGA X (Kumar et al. 
2018) using the Kimura two-parameter (K2P) model. The 
best-fit nucleotide substitution model was identified by the 
Akaike Information Criterion (AIC) in the jModelTest v.2 
(Posada 2008). AIC was TIM2+I+G for COI, GTR+I+G 
for Cyt b, GTR+I+G for the two mtDNA regions and 
TIM2+I+G for the mitochondrial and nuclear markers. 
Maximum Likelihood (ML) analysis was performed using 
RAXxMLy.7.2.6 (Stamatakis 2014) with 1000 replicate 
bootstrap. Maximum Likelihood (ML) analysis was 
performed using RAxMLv.7.2.6 (Stamatakis 2014) with 
1000 replicate bootstrap. Bayesian Inference (BI) was 
run in MrBayes v.3.2.7 (Ronquist et al. 2012) with two 
independent Markov Chain Monte Carlo for 50,000,000 
generations. The first 25% of trees were discarded as part 
of the burn-in phase. The remaining trees were combined 
to generate a 50% majority rule consensus tree. The 
posterior probability and bootstrap values were calculated 
from the BI and ML for each clade and values above 0.95 
and 75 were deemed as significant supports, respectively. 
A haplotype network was constructed using the median- 
joining method in the PopART v.1.7 package (Leigh and 
Bryant 2015). The parsimony informative and variable 
sites, as Well as haplotype diversity were calculated using 
DnaSP v.5.10 to test patterns of population structure in 
D. pictum and D. labiosum (Librado and Rozas 2009). 
Initially, genetic divergence between the putative 
species of Diagramma was estimated and the analysis 
of molecular variance (AMOVA) as ®,, values was 
performed in GenAIEx v.6.5 (Peakall and Smouse 2006) 
based on pairwise comparisons. Demographic changes in 
the effective population size were investigated for both 
Species using two neutrality tests (Fu’s Fs and Tajima’s 
D) and the Ramos-Onsins and Rozas test R2 in DnaSP 
v.5.10 (Librado and Rozas 2009). The time elapsed 
since most recent population expansion was calculated 
herein for Diagramma pictum and D. labiosum (Rogers 
and Harpending 1992) using the t = 2ut value, where t 
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is the time elapsed since the start of expansion in years, 
u is the mutation rate per generation for the sequence 
(the number of base pairs per divergence rate (1% or 2% 
Mya). generation time (in years) and t is the time elapsed 
since the beginning of expansion in years. As reported by 
Grandcourt et al. (2011), the generation time is 2.9 years 
for Diagramma pictum. The mutation rate reported for the 
mtDNA (COI) in previous reef fish studies ranged from 
1% to 2% per million years (Bowen et al. 2001). 


Divergence time estimation 


A calibrated time-tree of Plectorhinchinae was estimat- 
ed from the combined data set (COI, Cyt 6 and RAGI: 
3157 total bp) in BEAST v.2.6.2 (Bouckaert et al. 2019). 
Three independent MCMC chains were run for 20 million 
generations each with sampling every 1,000 generations. 
Runs were implemented using an uncorrelated lognormal 
relaxed clock model to calculate rate heterogeneity and 
the Yule prior method for the tree prior. Convergence of 
runs and the effective sample size (> 200) were checked 
in Tracer v.1.6 (Rambaut 2016). The three runs were 
combined by LogCombiner v.2.6.2 (Bouckaert et al. 
2019). The first 25% of trees were discarded in the Tree- 
Annotator v.2.6.2 (Bouckaert et al. 2019). The consensus 
tree with divergence time estimates was visualised using 
FigTree v.1.4.4 (Rambaut, 2016). Due to a lack of fossil 
record for Haemulidae, two calibration points were used 
following Near et al. (2013) and Tavera et al. (2018): (1) 
the Haemulidae-Lutjanidae split in early Eocene (mean 
52.5, 95% HPD interval 57.2-46.8) and (2) the Para- 
pristipoma-Plectorhinchus split in the late Eocene (mean 
34.5 Ma, 95% HPD interval 42.6—27 Ma). 


Ancestral area reconstruction 


The ancestral area in Plectorhinchinae was inferred using 
the BioGeoBEARS analysis in RASP v.4 (Yu et al. 2020). 
Geographic distribution of each species was derived from 
literature (Fricke et al. 2022) and sampling of our rang- 
es (Suppl. material 1). Four biogeographic regions were 
defined, based on present-day distribution and genetic 
diversification reported in previous studies (Ma et al. 
2016; Stiller et al. 2022): Western Atlantic (WA); Eastern 
Atlantic (EA); Western Indian Ocean (WIO); Central In- 
do-Pacific (CIP). The input trees of BEAST analysis were 
used to infer the ancestral area probability. Six different 
biogeographic models were evaluated to select the best- 
fit using the AIC cumulative weight (AICc_wt). 


Molecular species delimitation 


Molecular species delimitation was used from four dif- 
ferent methods, including; (A): the Automatic Barcode 
Gap Discovery (ABGD) (Puillandre et al. 2012) uses the 
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pairwise genetic distance to detect a barcode gap while as- 
suring that intra and interspecific genetic distances do not 
overlap and are executed in the website (http://wwwabi. 
snv.jussieu.fr/public/abgd) with prior intraspecific di- 
vergence (0.001 to 0.1), steps = 10 and a relative gap 
width of 1 using the Kimura 2-parameter model (Kimura 
(K80)); (B): Assemble Species by Automatic Partitioning 
(ASAP) (Puillandre et al. 2021), based on the distance 
threshold and a hierarchical clustering algorithm to cre- 
ate a list of best partitions. The partition score combines 
two metrics: the barcode gap width and the probability 
of panmixia (p-value), in which a lower score indicates 
a better partition. The method that was done using the 
website https://bioinfo.mnhn.fr/abi/public/asap with the 
Kimura (K80) model; (C): Bayesian Poisson Tree Pro- 
cess (bPTP) (Zhang et al. 2013), a method that estimates 
the number of mutations between two branching events 
using the branch length information of the input phyloge- 
ny. The bPTP approach was adopted on its website (http:// 
species. hits.org/ptp/) using a non-ultrametric tree of ML 
inference as the input, 500,000 MCMC generations and 
10% burn-in; (D): Generalised Mixed Yule Coalescent 
(GMYC) (Pons et al. 2006) method was run on its web- 
site (http://species.h-its.org/gmyc), using the ultrametric 
tree produced by BEAST v.2.6.2 (Bouckaert et al. 2019). 


Results 


MOTU delimitation and genetic diversity 


We assembled a total of 217 COI sequences (34 
newly developed plus 183 archived; 617 bp) from 29 
nominal Plectorhinchinae species. D. punctatum and 
Plectorhinchus pictus were sequenced for the first time. 
Novel COI barcodes were provided for Plectorhinchus 
flavomaculatus, P. gaterinus and P. gibbosus for the first 
time from the Northern Indian Ocean (NIO). Of 29 nominal 
species studied, the coalescent-based delimitations (bPTP, 
GMYC) yielded 27 MOTUs, whereas the distance-based 
methods (ABGD and ASAP) yielded 26 MOTUs (Fig. 1). 
ABGD and ASAP identified Plectorhinchus lessonii, and 
P. diagrammus and P. albovittatus and P. caeruleonothus 
as a MOTU. bPTP and GMYC methods considered 
Plectorhinchus lessonii and P. diagrammus as one MOTU. 
Additionally, haplotype networks, based on the COI gene, 
identified two main distinct haplogroups (Fig. 2A). The 
haplogroup I was composed of Diagramma pictum from 
the NIO to WP, which was separated by 15 mutation steps 
from haplogroups II, represented by D. punctatum from 
the Red Sea. In the first haplogroups, samples of D. pictum 
were mixed with haplotypes of D. Jabiosum from 
Australia. As suggested by the genetic distance method 
that uses the K2P model, the highest interspecific genetic 
distance was found between Plectorhinchus macrolepis 
and P. schotaf (21.3%). In comparison, the lowest distance 
was observed between D. pictum and D. labiosum (0.53%) 
and P. lessonii and P. diagrammus (1.2%) (Table 1). The 
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Table 1. Genetic K2P distances (%) of the Diagramma and 
Plectorhinchus species. 


Species Maximum intra- Nearest- Distance 
group distance Neighbour (NN) to NN 
Diagramma pictum 0.17 D. labiosum 0.53 
D. punctatum 0.00 P. centurio 3.45 
Plectorhinchus 0.00 P. caeruleonothus 25a 
albovittatus 
P. caeruleonothus 0.19 P. albovittatus Ask 
P. centurio 0.35 D. punctatum 3.45 
P. chaetodonoides 0.43 P. centurio 8.75 
P. chrysotaenia 0.00 P. playfairi 16.30 
P. chubbi 0.08 P. sordidus B52 
P. cinctus 0.64 P. gibbosus 13.44 
P. diagrammus 0.00 P. lessonii 1.20 
P. flavomaculatus 0.31 P. makranensis 12.62 
P. gaterinus 0.28 P. caeruleonothus 12.18 
P. gibbosus 0.69 P. cinctus 13.44 
P. lessonii 0.08 P. diagrammus 1.20 
P. lineatus 0.06 P. polytaenia 7.17 
P. macrolepis 0.00 P. plagiodesmus 13.84 
P. macrospilus 0.00 P. lessonii 4.10 
P. makranensis 0.01 P. schotaf owas: 
P. orientalis - P. vittatus 7.97 
P. pictus 0.00 P. cinctus 13:55 
P. picus 0.08 D. punctatum 8.38 
P. plagiodesmus 0.00 P. pictus 13.24 
P. playfairi 0.16 P. chubbi 4.61 
P. polytaenia 0.00 P. lineatus Pairs 
P. schotaf 0.13 P. makranensis 523 
P. sordidus 0.10 P. chubbi 3.52 
P. unicolor 0.23 P. chubbi 5.80 
P. vittatus - P. orientalis 7.97 


highest intraspecific genetic variation was observed in 
P. gibbosus (0.69%), while values close to 0 to 0.09% were 
detected in most of the Diagramma and Plectorhinchus 
species (Table 1). Based on all delimitation analyses, 
D. labiosum specimens from Queensland, Australia (type 
locality) are genetically analogous to D. pictum and there 
is no evidence of their genetic differentiation. Cyt b 
sequence of a single D. Jabiosum specimen (DQ784750), 
deposited at the GenBank of Australia with D. pictum 
specimens in the phylogenetic tree, was identified as a 
monophyletic lineage (genetic divergence: 0.51%) (data 
are not shown). Seventeen haplotypes were detected based 
on COI sequences between D. labiosum and D. pictum 
(Table 2). Seventeen haplotypes were detected, based 
on COI sequences between D. labiosum and D. pictum 
(Table 2). The highest haplotype and nucleotide diversity 
indices were found in the samples of D. /abiosum 
(Australia) (Hd = 0.600, 2 = 0.0021), while specimens of 
D. pictum (Hd =0.560: NIO, 7=0.0018: WP) demonstrated 
the lowest values (Hd = 0.579, 7 = 0.001) (Table 2), where 
the diversity indicators were not influenced by sample 
size reduction. OST value revealed a non-significant 
level between Australia (D. labiosum) and IWP, except 
Australia (D. pictum) (®,. = —1.783) (Table 2). Test of 
neutrality found negative and non-significant Tajima’s 
D values, while Fu’s Fs (-8.801) and R2 (0.038) results 
were significant (Table 2). These significant values (Fu’s 
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Figure 1. Bayesian Inference tree, based on the COI showing the results of species delimitation algorithms. The vertical bars in- 
dicate species delimitation analyses. The black circles indicate nodes with supports (ML bootstrap BP > 75% and BI probability 
PP => 0.95), orange circles (BP < 75% and BI PP = 0.95) and blue circles (PP < 0.95 and ML BP < 75%). 


Fs, R2), in combination with the estimates of t for both 
Species, indicate a shallow population history with recent 
expansion (Table 2). The unimodal mismatch distribution 
for the pooled sample of D. pictum supported a recent 
historical expansion (Fig. 3). 


Phylogenetic relationships 


We analysed 254 COI DNA sequences (617 bp) which 
show 325 variable sites; the Cyt b gene consisting of 1117 
bp showed 295 variable positions; and Rag1 consisting 
of 1440 bp showed 184 variable sites. Neither gaps nor 
stop codons were observed in the final alignment. ML 
and BI analyses with the combined dataset (mtDNA and 
nuclear) revealed that the phylogenetic relationships of 
sweetlips species had the highest resolution (Fig. 4). The 
ML and BI phylogenetic trees, based on two mtDNA (Cyt 
b + COJ) and one nuclear (Ragl) gene, had a congruent 


topology, with a slight difference in their support values. 
The phylogenetic relationships observed in the analysis 
of the concatenated data (Cyt b + COI + Rag1) revealed 
that Parapristipoma, as a basal position, is the sister clade 
of the Plectorhinchus species with high bootstrap and 
posterior probabilities (Fig. 4). Within the Plectorhinchus 
clade, two well-supported subclades were formed: (1) The 
first clusters the most WIO species (subclade I): Plecto- 
rhinchus schotaf, P. makranensis, P. playfairi, P. sordi- 
dus, P. chubbi, P. unicolor, P. plagiodesmus, P. gibbosus, 
P. macrolepis, P. chrysotaenia, P. cinctus and P. pictus; 
(2) the second clusters the most CIP species: P. gaterinus, 
P. chaetodonoides, P. picus, P. albovittatus, P. lessonii, 
P. lineata, P. macrospilus, P. diagrammus, P. polytaenia, 
P. vittatus, P. orientalis, P. caeruleonothus along with 
D. pictum and D. punctatum indicating that the Diagram- 
ma species is nested within Plectorhinchus and should be 
recognised as Plectorhinchus pictum and Plectorhinchus 
punctatum (subclade IT) (PP = 1; BS = 100%). 
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Figure 2. Median-joining network obtained of 617 bp of mitochondrial COI for Diagramma species (A). Map depicting the geo- 
graphical distribution of 17 haplotypes (B). Circle sizes are relative to haplotype frequency. Numbers between haplotypes represent 
mutational steps between haplotypes. Colour circles indicate species of D. pictum (red: Persian Gulf to Bangladesh, green: China 
and Indonesia) and D. /abiosum (orange: Australia), D. punctatum (black: Red Sea). 
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Pairwise differences 
Figure 3. Mismatch distribution inferred from partial mtDNA 
COI for the pooled dataset of D. pictum. 


Divergence times and ancestral area 
reconstructions 


The Haemulidae diverged from Lutjanidae ~ 49 Ma (ear- 
ly to middle Eocene, 95% HPD: 43.6—55.2 Ma) and then 
Plectorhinchinae divided into two main clades (A-—B) 
and two subclades (I-II) (Fig. 5). The first split created 
Clade A containing Parapristipoma at ~ 35 Ma (late Eo- 
cene, 95% HPD: 32.4—37.1 Ma). Clade B diverged at ~ 
21 Ma (early Miocene, 95% HPD: 18.4—23.5 Ma) into 
two subclades. Within each subclade, there were two pe- 
riods of maximum diversification from the Miocene to 
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Table 2. Molecular diversity indices for D. pictum and D. labio- 
sum based on mitochondrial DNA (COI) sequence data. Numbers 
in bold denotes significance (P< 0.05) and for Fu’s Fs (P< 0.02). 


Statistics D. pictum D. labiosum 
Japan to Iran to Australia All 
Indonesia Bangladesh (WP) samples 
(WP) (NIO) 
Number of individuals (N) N = 44 N = 26 N=/7 N = 76 
Expansion time (years) _ - - 17,647- 
35,294 
Number of haplotypes (Hn) 8 6 3 17 
Haplotype diversity 0.579 + 0.560 + 0.600+ 0.624+ 
(Hd + SD) 0.075 0.073 0.215 0.051 
Nucleotide diversity 0.0018+ 0.0021+ 0.0021+ 0.0024 + 
(7 + SD) 0.0011 0.0016 0.0017 0.0018 
Tajima’s D -1.218 -1.636 -1.233 -1.783 
Fu’s FS -3.739 -1.067 0.189 -8.801 
R2 0.058 0.060 0.062 0.038 
O., Australia 0.084 0.072 = 
lran to 0.049 - - 
Bangladesh 


Pliocene epochs. The best-fit biogeographic models re- 
covered for the ancestral range reconstruction was the 
BAYAREALIKE (AICc_wt: 0.71). As suggested by the 
ancestral area reconstruction using the combined dataset 
(Cyt b + COI + Ragl1), 17 dispersal and eight vicariance 
events occurred during the evolution of the studied Plec- 
torhinchinae. The deepest node of early Plectorhinchinae 
constructed probably has its origin in the [TWP with the 
dispersal and vicariance events occurring from the late 
Eocene to early Oligocene epoch (Fig. 6). However, more 
comprehensive taxonomic sampling of the Parapristipo- 
ma species can help find a definitive point for the origin 
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Figure 4. Phylogeny recovered by the Bayesian Inference (BI) and Maximum Likelihood (ML) analyses for the Cyt b, COI and 
RAGI dataset. The black circles indicate nodes with supports (ML bootstrap BP = 75% and BI probability PP = 0.95), orange circles 
(BP < 75% and BI PP = 0.95) and blue circles (PP < 0.95 and ML BP < 75%). 


of the subfamily Plectorhinchinae. The ancestors of the 
earliest Plectorhinchus lineages were supposed to be of 
WIO origin with vicariance and dispersal events taking 
place in the Miocene. Within subclade I, there seems to 
be a burst of species diversification in the WIO at 20 Mya. 
Within subclade II, the vicariance event affected specia- 
tion in the CIP as the ancestral area at 15 Mya. 


Discussion 


This study represents the first attempt to delimit the 
species boundaries and evolutionary and biogeographic 
history of marine fish species in Plectorhinchinae. The 
present study corroborates the findings of the previous 
molecular studies that show paraphyly of the Plecto- 
rhinchus without the inclusion of Diagramma (San- 
ciangco et al. 2011; Tavera et al. 2018) and contrast 
with morphological data (Konchina 1976; Carpenter 
and Niem 2001; Nelson et al. 2016) that supported this 
genus monophyly. 


Molecular phylogeny of Plectorhinchus 


According to the molecular data, Plectorhinchus is divid- 
ed into two subclades (Fig. 4), which display high concor- 
dance with the division proposed by Tavera et al. (2018). It 
can reinforce the phylogenetic value of the geographic dis- 
tribution and body-colour pattern in distinguishing P/ecto- 
rhinchus lineages. The first group (subclade I) mostly con- 
sists of species with uniformly dark colours (Plectorhinchus 
schotaf, P. makranensis, P. playfairi, P. sordidus, P. chub- 
bi, P. unicolor, P. plagiodesmus, P. gibbosus, P. macrolepis 
and P. chrysotaenia) or black pots on the upper sides of the 
body (P. cinctus and P. pictus), along with more profound 
body depth (Randall 1995, 1997; McKay 2001), which are 
those typically distributed in the WIO (Fig. 4) (Fricke et al. 
2022). The second group (subclade II) was comprised of 
species with spots on the whole body (P. gaterinus, P. cha- 
etodonoides and P. picus) or horizontal stripes alongside 
(P. albovittatus, P. lessonii, P. lineata, P. macrospilus, 
P. diagrammus, P. polytaenia, P. vittatus and P. orienta- 
lis), except for P. caeruleonothus and Diagramma species 
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Figure 5. Time-calibrated phylogeny, based on combined mitochondrial (COI, Cyt 5) and nuclear (RAG1) genes, inferred from the 


BEAST analysis of the Plectorhinchinae. Numbers in the horizontal purple bars represent the 95% Highest Posterior Density inter- 


vals for each estimated event. The black circles (pp = 0.95) and orange circles (pp < 0.95). A: Parapristipoma, B: Plectorhinchus. 


(Fig. 4) (Randall 1995, 1997; Johnson et al. 2001; McKay 
2001). The species includes representatives with shorter 
bodies from the CIP (Fricke et al. 2022). The colour varia- 
tion observed in subclade II compared to subclade I can be 
attributed to the substrate type (Randall et al. 1997; Mck- 
ay 2001; Johnson et al. 2001). The colour variation of the 
members of subclade II inhabiting the coral reef may be 
linked to social communication, mimicry and camouflage 
defence mechanisms to avoid predation adjusted to their 
habitat (Randall 1998; McMillan et al. 1999). 


Systematics and species boundaries of 
Plectorhinchus 


In our analyses, specimens from 29 nominal species were 
delimited between 26 and 27 consensus MOTUs. In gen- 
eral, low interspecific genetic divergences were observed 
in some of the MOTUs within Plectorhinchus, which 
focuses on the barcoding analysis in closely-related spe- 
cies. This situation was formerly reported in haemulid 
within Anisotremus (Bernardi et al. 2008) and Haemu- 
Jon (Rocha et al. 2008). In the meanwhile, the fact that 
P. lessonii and P. diagrammus, along with P. albovittatus 
and P. caeruleonothus, grouped together in all the delim- 
itation analyses, exhibits their very recent divergence (< 
1 Ma) compared to other Plectorhinchus members. The 
morphological variability, including minor differences in 
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colouration between D. pictum and D. labiosum deter- 
mines taxonomic classification into subspecies or species 
status (Allen 1997; Johnson et al. 2001; Bogorodsky et 
al. 2014; Bray 2017; Fricke et al. 2022). Our molecular 
analyses did not yield any genetic evidence to confirm D. 
labiosum as a distinct species (Figs 1, 2, Table 1). A rela- 
tively low level of genetic divergence between D. pictum 
and D. labiosum found in both markers (COI and Cyt 5) 
may indicate a recent speciation event (Table 1). These 
distances are considerably shorter than those reported in 
haemulid species (Johnson and Wilmer 2015; Damadi et 
al. 2020) and other related species of marine fishes, which 
is at least 2% between closely-related species (Ward et al. 
2005; Hsu et al. 2007). Demographic tests (Fu’s, R2 and 
mismatch distribution) indicate that D. Jabiosum and D. 
pictum share a common ancestor and expansion events 
over the last few thousand years coincided with the Ho- 
locene glaciations that restricted larval dispersal and 
ocean circulation deflection (Yokoyama et al. 2001; Wil- 
son 2013; Rosser et al. 2020). These features might have 
reduced coral reef habitats and colonisation from north- 
west Australia, especially in the reef fish (Ovenden et al. 
2002; Salini et al. 2006; Sulaiman and Ovenden 2010). 
The incongruence between colour variation and the lev- 
el of genetic divergence is not rare in fish and could be 
attributed to some other reasons: (1) phenotypic plastici- 
ty rather than reproductive isolation, (2) incomplete lin- 
eage sorting and (3) introgression through hybridisation 
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Figure 6. Ancestral area reconstruction obtained using BioGeoBEARS for Plectorhinchinae. Biogeographic areas are shown in the 


circles and the colours in the squares indicate the current species distribution. 


(Bowen et al. 2006; Schultz et al. 2007; DiBattista et al. 
2012; DiBattista et al. 2015; Bernal et al. 2017). There are 
several lines of evidence that favour the incomplete lin- 
eage sorting over other scenarios; however, there is still 
a paucity of evidence and further analyses using genomic 
data are warranted. Firstly, D. Jabiosum and D. pictum 
showed extensive mtDNA sharing and coalescence of 
haplotypes, which could be interpreted as evidence of a 
very recent divergence (Figs 1, 2). Secondly, their distri- 
bution is allopatric (D. /abiosum: Australia and D. pic- 
tum: IWP, except Australia) (Fig. 2B) (Fricke et al. 2022). 
Hybridisation events and lineage sharing usually occur 
in sympatric and parapatric areas rather than in allopatric 
areas due to the greater opportunities for hybridisation in 
the secondary contact zones (Hobbs et al. 2009; DiBattis- 
ta et al. 2016). In contrast, incomplete lineage may take 
place in both sympatric and allopatric areas (DiBattista et 
al. 2012; DiBattista et al. 2015). Thirdly, regional colour 
morphs usually represent minor molecular divergences 
due to recent genetic divergence, which may not be re- 
lated to ecological and behavioural adaptations (Taylor 
and Hellberg 2003; Drew et al. 2008; Drew et al. 2010). 
Fourthly, the Diagramma species are distinguished only 
based on colour patterns (Johnson et al. 2001; Heems- 
tra and Heemstra 2004; Bogorodsky et al. 2014) which 
evolve more quickly than the other morphological traits 
and genetic characters (Bowen et al. 2006; Schultz et 
al. 2007). Colouration could be an important criterion 


to distinguish closely-related species in haemulid fishes 
(Johnson and Wilmer 2015; Damadi et al. 2021). Fifthly, 
the most recent divergence occurred between D. pictum 
and D. labiosum around 28—56 million years ago (Table 
1), which is more recent than the complete lineage sorting 
time. A divergence time of 2.9 Mya seems sufficient for 
these lineages for complete sorting (1 generation time) 
(Grandcourt et al. 2011). 


Biogeography and temporal diversification 


Our biogeographic findings exhibit a complex history 
highlighting the role of both biological and geographical 
factors in the diversification and evolution within the Hae- 
mulids. The first diversification, 1.e. the split of Haemulids 
from Lutjanids as a sister lineage in the Eocene, could be 
related to the ice sheet formed in the Antarctic and reduced 
atmospheric CO, level (Zachos et al. 2001; Coxall et al. 
2005; Lear et al. 2008). These climatic and oceanographic 
changes resulted in a sea-level drop, upwelling and extinc- 
tion of warmer climate fauna (Zachos et al. 2001; Liver- 
more et al. 2007; Ma et al. 2016). This extinction formed 
vacant niches and the emerging Indo-Australian Archipel- 
ago probably gave rise to new reef habitats and high spe- 
cies-richness in the cooler oceans (Cowman and Bellwood 
2011; Ma et al. 2016). This divergence time corresponds 
to the one shown by Near et al. (2012) for the same his- 
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torical event in Haemulids. The historical biogeography 
suggests that the common ancestor of Plectorhinchinae 
might have its origin in the [WP and subsequently spread 
and diversified in the western Indian Ocean and the Cen- 
tral Indo-Pacific. This scenario could be firmly demon- 
strated by complete biogeographic sampling of some taxa 
and extra genetic/genomic analyses. For Plectorhinchinae, 
divergence was first observed in East Atlantic (EA) spe- 
cies (Parapristipoma) with the fewest species from the 
Plectorhinchus \ineage, approximately 33.8 Mya, which 
coincided with the onset of the terminal Tethyan event 
(TTE) that disconnected West Tethyan and Arabian ma- 
rine biodiversity hotspots (Adams et al.1983; Harzhauser 
et al. 2007). The diminishing connectivity between two re- 
gions might have set the scene for vicariant splits, which 
appear to be linked to the early diversification of this clade. 
This divergence time is also congruent with the EA and 
IWP divergence in other species of reef fish, particularly 
the Groupers and Boxfishes (Santini et al. 2013; Ma et al. 
2016). Within the Plectorhinchus TWP Clade, two factors, 
the Indo-Australian Archipelago (IAA) uplift and major 
global cooling, were associated with the diversification 
of sweetlips subclades. These factors were inferred from 
other reef fish in the Miocene period (Zachos et al. 2001; 
Bellwood and Wainwright 2002; Shevenell et al. 2004; 
Cowman and Bellwood 2011; Ma et al. 2016). An esti- 
mate of divergence times of each Plectorhinchus subclade 
suggests further diversification occurred during the Mio- 
cene that slow indicated by the long branches, followed 
by slower diversification in Pliocene to Pleistocene epochs 
(Fig. 5). The IAA uplift during the early to mid-Miocene 
in subclade II (CIP Plectorhinchus) created new islands, 
coinciding with the lowering of sea level, leading to novel 
shallow reef habitat discontinuity. During Plio-Pleistocene, 
Diagramma species diverged from each other around 4.1 
to 2.3 Mya (Fig. 5) and D. punctatum, an endemic spe- 
cies, formed in the Red Sea, coinciding probably with high 
temperature and salinity, sea-level changes during glacial 
cycles and the narrow Strait of Bab al Mandab. These fac- 
tors led to a decrease in unsuitable habitat and the water 
flow through this strait (Rihm and Henke 1998; Riegl and 
Purkis 2012; Stevens et al. 2014; Torquato et al. 2019) that 
has been documented in marine fish (Froukh and Kochzius 
2007; DiBattista et al. 2013; Bogorodsky et al. 2017). The 
major global cooling could lay the ground for emergence 
of new environmental niches likely driven by adaptive ra- 
diations in subclade I (WIO) during the Miocene (Zachos 
et al. 2001; Shevenell et al. 2004). 


Conclusions 


The findings of the present study on _ sweetlips 
(Plectorhynchinae) not only shed light on the interrelation- 
ships and evolutionary history, but also opened new avenues 
for research. The results of this study provide a framework 
to shed light on the probable origin and factors involved 
in the diversification of the [WP Plectorhynchinae. The 
contemporary species have a common ancestor that comes 
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from the present day [WP and the current geographical dis- 
tribution appears to be related to vicariance and dispersal 
events. It indicates that Diagramma species were divided 
4.1 to 2.3 Mya ago into ancestral Diagramma punctatum 
in the Indian Ocean, directly followed by the WP and re- 
cently isolated species. These findings raise the possibility 
that recent genetic divergence have given rise to colour 
morphotypes as emerging species on coral reefs. 
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